Method for estimating nondirectional wave spectrum from sea echoes of multiple high radar frequencies

ABSTRACT

The disclosure provides a method for estimating the nondirectional wave spectrum from the sea echoes of multiple HF radar frequencies. The method includes: dividing the radar detection area into a plurality of fan-shaped units at an equal range interval and angle interval according to the distance resolution and the angular resolution of an HF radar; obtaining the Doppler spectrum from the sea echo of a single radar frequency at a fan-shaped unit by performing the first fast Fourier transform (FFT) in distance dimension, the second FFT in Doppler frequency dimension and the digital beamforming; extracting the positive first-order peak and the negative first-order peak from the aforementioned Doppler spectrum by the peak-searching method; and selecting the stronger first-order peak σR(1)(ω) ; dividing the second-order spectrum on the stronger first-order peak side into an inner second-order spectrum and an outer second-order spectrum.

CROSS-REFERENCE TO RELAYED APPLICATIONS

Pursuant to 35 U.S.C.§ 119 and the Paris Convention Treaty, thisapplication claims foreign priority to Chinese Patent Application No.202110352889.X filed Apr. 1, 2021, the contents of which, including anyintervening amendments thereto, are incorporated herein by reference.Inquiries from the public to applicants or assignees concerning thisdocument or the related applications should be directed to: MatthiasScholl P. C., Attn.: Dr. Matthias Scholl Esq., 245 First Street, 18thFloor, Cambridge, Mass. 02142.

BACKGROUND

The disclosure relates to a method for remote sensing of the oceansurface, and more particularly to a method for estimating thenondirectional wave spectrum from the sea echoes of multiple highfrequency (HF) radar frequencies.

In recent years, HF radars are widely used in marine environmentalmonitoring, which has the characteristics of extensive detectioncoverage, all-weather, all-day and high measurement accuracy. HF radarstransmit the vertical polarization electromagnetic waves that propagatealong the ocean surface with little attenuation. The electromagneticwaves interact with the ocean surface to generate Bragg scattering, andthen the Doppler spectrum containing the information of the sea statecan be obtained. Wave parameters can be retrieved from the Dopplerspectrum, including the ocean wave spectrum, the wave height and thewave period, etc. Most existing techniques estimate the wave informationby using single-frequency HF radars. Under a high sea state, thesecond-order part of the Doppler spectrum is prone to be saturated anddifficult to be separated from the first-order peak for a high radarfrequency. Nevertheless, under a low sea state, the energy of thesecond-order part is usually low and submerged in the noise floor for alow radar frequency. Therefore, the extraction of wave parameters byusing a single-frequency HF radar is not robust enough for themulti-scale wave detection under different sea states, because themeasurable range of wave height from the single-frequency data islimited. To broaden the measurable range of wave height and improve therobustness of the measurement under different sea states, it isnecessary to directly estimate the nondirectional wave spectrum from thesea echoes of multiple HF radar frequencies.

SUMMARY

Accordingly, it is an object of the disclosure to provide a method forestimating the nondirectional wave spectrum from the sea echoes ofmultiple HF radar frequencies.

It is another object of the disclosure to provide such a method to solvethe problem of limited measurable range of wave height observed at afixed radar frequency by combining echo signals from multiple HF radarfrequencies.

It is a further object of the disclosure to provide such a method toimprove the accuracy of measuring the nondirectional wave spectrum andmeet the requirements of robustly extracting wave parameters undervarious sea states.

The disclosure relates to a method for estimating the nondirectionalwave spectrum from the sea echoes of multiple HF radar frequencies. Themethod is carried out in an apparatus which comprises a multi-frequencyHF radar and a computer processor. The multi-frequency HF radarcomprises a transmitter for transmitting electromagnetic wave signalswith multiple frequencies between 3-30 MHz in the HF band to the oceansurface; a receiver for receiving sea echoes of multiple radarfrequencies originated from the electromagnetic wave signals modulatedby ocean surface movements; and a signal processor unit for obtainingeach Doppler spectrum from the sea echo of each radar frequency at afan-shaped unit by performing a first fast Fourier transform (FFT) indistance dimension, a second FFT in Doppler frequency dimension and adigital beamforming. The computer processor provides code segments forestimating the nondirectional wave spectrum.

The method for estimating the nondirectional wave spectrum from the seaechoes of multiple HF radar frequencies comprises:

(a) dividing the radar detection area into a plurality of fan-shapedunits at an equal range interval and angle interval according to adistance resolution and an angular resolution of a multi-frequency HFradar, wherein the multi-frequency HF radar is capable of simultaneouslyoperating at more than one frequency in the HF band;

(b) obtaining a Doppler spectrum from a sea echo of a single radarfrequency at a fan-shaped unit by performing a first fast Fouriertransform (FFT) in distance dimension, a second FFT in Doppler frequencydimension and a digital beamforming in a signal processor unit,extracting a positive first-order peak and a negative first-order peakfrom the Doppler spectrum by a peak-searching method and selecting astronger first-order peak σ_(R) ⁽¹⁾(ω) from the positive and negativefirst-order peaks in a computer processor;

(c) dividing a second-order spectrum on a stronger first-order peak sideinto an inner second-order spectrum and an outer second-order spectrumand separating the outer second-order spectrum on the strongerfirst-order peak side from the Doppler spectrum in the computerprocessor according to a Doppler frequency range of the outersecond-order spectrum;

(d) calculating R_(f)(ω) which is a function of the outer second-orderspectrum divided by a first-order peak energy for the stronger side ofthe Doppler spectrum of the single radar frequency at the fan-shapedunit in the computer processor;

(e) calculating A_(f) which is defined as a single-frequency coefficientmatrix of the nondirectional wave spectrum by linearizing the outersecond-order spectrum on the side of the stronger first-order peak inthe computer processor;

(f) repeating the above steps (b) to (e), combining the sea echoes ofdifferent radar frequencies to construct a matrix R which is a functionof the outer second-order spectrum divided by the first-order peakenergy for the stronger side of the Doppler spectrum of multiple radarfrequencies and merging coefficient matrices of the nondirectional wavespectrum from multiple radar frequencies into a matrix A in the computerprocessor; and

(g) calculating a pseudo-inverse A⁺ of the matrix A by a singular valuedecomposition and estimating the nondirectional wave spectrum at thefan-shaped unit from the matrix R and the pseudo-inverse A⁺ in thecomputer processor.

In a class of this embodiment, the Doppler spectrum at the fan-shapedunit stated in the step (b) is defined as: σ(ω), where co represents aDoppler frequency generated by the motion of ocean waves to themulti-frequency HF radar; and σ(ω) represents a wave energy distributionat the different value of ω.

The positive first-order peak and negative first-order peak areextracted from the Doppler spectrum at the fan-shaped unit using thepeak searching method. The first-order peaks are defined as two peaks inthe Doppler spectrum which are roughly symmetrically distributed on bothsides of zero frequency. The first-order peaks are generated by theBragg scattering of the waves of half the radar wavelength that areeither advancing directly towards the multi-frequency HF radar orreceding directly from the multi-frequency HF radar.

The step (b) comprises the steps of:

searching for a point of a maximum amplitude in the Doppler frequencyrange [0.6ω_(B),1.4ω_(B)] of the Doppler spectrum at the fan-shaped unitas the peak of the positive first-order peak, and recording the Dopplerfrequency of the peak of the positive first-order peak as ω_(P+),wherein ω_(B)=√{square root over (2gk₀)} is a Bragg frequency, wherein

$k_{0} = \frac{2\pi f}{c}$

is a radar wavenumber; c is a speed of light; and f is a radarfrequency;

searching for a local minimum point inside the peak of the positivefirst-order peak where the Doppler frequency meets the inequationω_(P+)−0.2ω_(B)≤ω<ω_(P+), and denoting the Doppler frequency of thelocal minimum point inside the peak of the positive first-order peak asω_(L+);

searching for the local minimum point outside the peak of the positivefirst-order, peak where the Doppler frequency satisfies the inequationω_(P+)<ω≤0.2ω_(B)+ω_(P+), and recording the Doppler frequency of thelocal minimum point outside the peak of the positive first-order peak asω_(R+);

intercepting the Doppler spectrum at the fan-shaped unit with Dopplerfrequency [ω_(L+),ω_(R+)] as the positive first-order peak;

searching for the point of the maximum amplitude in the Dopplerfrequency range [−1.4ω_(B), −0.6ω_(B)] of the Doppler spectrum at thefan-shaped unit as the peak of the negative first-order peak, andrecording the Doppler frequency of the peak of the negative first-orderpeak as ω_(P−);

searching for the local minimum point inside the peak of the negativefirst-order peak where the Doppler frequency meets the inequationω_(P−)<ω≤ω_(P−)0.2ω_(B), and denoting the Doppler frequency of the localminimum point inside the peak of the negative first-order peak asω_(L−);

searching for the local minimum point outside the peak of the negativefirst-order peak where the Doppler frequency satisfies the inequationω_(P−)−0.2ω_(B)≤ω<ω_(P−), and recording the Doppler frequency of thelocal minimum point outside the peak of the negative first-order peak asω_(R−);

intercepting the Doppler spectrum at the fan-shaped unit with Dopplerfrequency [ω_(R−),ω_(L−)] as the negative first-order peak; and

comparing the amplitude of the peak of the positive first-order peakwith that of the negative first-order peak, and selecting thefirst-order peak with a larger amplitude of the peak point as thestronger first-order peak σ_(R) ⁽¹⁾(ω).

In a class of this embodiment, in the step (c), the second-orderspectrum originated from the second-order scattering of ocean waves andradar waves is a continuum with lower amplitude than the first-orderpeaks and distributed around the first-order peaks.

The step (c) comprises the steps of:

(1) dividing the second-order spectrum on the stronger first-order peakside into an inner second-order spectrum and an outer second-orderspectrum;

The Doppler frequency range of the outer second-order spectrum on thestronger first-order peak side is given as

$\left\{ \begin{matrix}{{\omega_{c +} < \omega \leq {1.4\omega_{B}}},{{when}{\sigma_{R}^{(1)}(\omega)}{is}{the}{positive}{first} - {order}{peak}}} \\{{{{- 1.4}\omega_{B}} \leq \omega_{c -}},{{when}{\sigma_{R}^{(1)}(\omega)}{is}{the}{negative}{first} - {order}{peak}}}\end{matrix} \right.$

The Doppler frequency range of the inner second-order spectrum on thestronger first-order peak side is given as

$\left\{ \begin{matrix}{{{0.6\omega_{B}} \leq \omega < \omega_{c +}},{{when}{\sigma_{R}^{(1)}(\omega)}{is}{the}{positive}{first} - {order}{peak}}} \\{{\omega_{c -} < \omega \leq {- 0.6\omega_{B}}},{{when}{\sigma_{R}^{(1)}(\omega)}{is}{the}{negative}{first} - {order}{peak}}}\end{matrix} \right.$

where ω is a Doppler frequency generated by the motion of ocean waves toa radar; ω_(B)=√{square root over (2gk₀)} is a Bragg frequency; k₀ is aradar wavenumber;

$\begin{matrix}{\omega_{c +} = \frac{\sum\limits_{\omega_{L +}}^{\omega_{R +}}{\omega \cdot {\sigma_{R}^{(1)}(\omega)}}}{\sum\limits_{\omega_{L +}}^{\omega_{R +}}{\omega \cdot {\sigma_{R}^{(1)}(\omega)}}}} \\

\end{matrix}$ isacentroidfrequencyofthepositivefirst − orderpeak; and$\omega_{c -} = \frac{\sum\limits_{\omega_{R -}}^{\omega_{L -}}{\omega \cdot {\sigma_{R}^{(1)}(\omega)}}}{\sum\limits_{\omega_{R -}}^{\omega_{L -}}{\omega \cdot {\sigma_{R}^{(1)}(\omega)}}}$isacentroidfrequencyofthenegativefirst − orderpeak.

(2) if the stronger first-order peak obtained in the step (b) is thepositive first-order peak, taking the centroid frequency ω_(c+) as astarting point of Doppler frequency, taking a cutoff frequency of1.4ω_(B) as an endpoint, extracting the Doppler spectrum at thefan-shaped unit σ(ω) in the frequency range of (ω_(c+),1.4 ω_(B)], andsetting the Doppler spectrum at the fan-shaped unit σ(ω) in thefrequency range of (ω_(c+), ω_(R+)] to zero to obtain the outersecond-order spectrum on the stronger first-order peak side σ_(R)⁽²⁾(ω);

if the stronger first-order peak obtained in the step (b) is thenegative first-order peak, taking the centroid frequency ω_(c−) as thestarting point of Doppler frequency, taking a cutoff frequency of−1.4ω_(B) as the endpoint, extracting the Doppler spectrum at thefan-shaped unit σ(ω) in the frequency range of [−1.4ω_(B),ω_(c−)), andsetting the Doppler spectrum at the fan-shaped unit σ(ω) in thefrequency range of [W_(R−), co) to zero to obtain the outer second-orderspectrum on the stronger first-order peak side σ_(R) ⁽²⁾(ω).

Preferably, in the step (d), the aforementioned function R_(f)(ω) isgiven by the equation:

$\left\{ \begin{matrix}{{{R_{f}(\omega)} = \frac{\sigma_{R}^{(2)}(\omega)}{\sum\limits_{\omega_{L +}}^{\omega_{R +}}{{\sigma_{R}^{(1)}(\omega)}\Delta\omega}}}\ ,\ {{when}\ {\sigma_{R}^{(1)}(\omega)}\ {{is}\ {the}\ positive\ first}}‐{order\ peak}} \\{{{R_{f}(\omega)} = \frac{\sigma_{R}^{(2)}(\omega)}{\sum\limits_{\omega_{R -}}^{\omega_{L -}}{{\sigma_{R}^{(1)}(\omega)}{\Delta\omega}}}}\ ,\ {{when}\ {\sigma_{R}^{(1)}(\omega)}\ {{is}\ {the}\ negative\ first}}‐{order\ peak}}\end{matrix} \right.$

where f in a subscript indicates a radar frequency; σ_(R) ⁽¹⁾(ω) is thestronger first-order peak; σ_(R) ⁽²⁾(ω) is the outer second-orderspectrum on the stronger first-order peak side; ω_(L+) is the Dopplerfrequency of the local minimum point inside the peak of the positivefirst-order peak; ω_(R+) is the Doppler frequency of the local minimumpoint outside the peak of the positive first-order peak; ω_(L−) is theDoppler frequency of the local minimum point inside the peak of thenegative first-order peak; ω_(R−) is the Doppler frequency of the localminimum point outside the negative first-order peak; and Δω is afrequency resolution of the Doppler spectrum σ(ω).

In a class of this embodiment, the step (e) comprises the steps of:

(1) linearizing the outer second-order spectrum on the side of thestronger first-order peak;

at the Doppler frequencies near the first-order peak that satisfy thecondition ω_(B)<|ω|≤1.4ω_(B), the direction of one of the two sets ofocean waves that generate the second-order scattering with the radarvector is approximately equal to that of the Bragg wave vector, so theouter second-order spectrum σ_(R) ⁽²⁾(ω) can be linearized as σ_(RL)⁽²⁾(ω):

${{\sigma_{RL}^{(2)}(\omega)} = {\sum\limits_{\theta = 0}^{2\pi}{2^{8}\pi k_{0}^{4}{❘\Gamma ❘}^{2}{S\left( {m\overset{\rightarrow}{k}} \right)}{S\left( {{- 2}m^{\prime}\overset{\rightarrow}{k_{0}}} \right)}\frac{\left( {2k_{0}} \right)^{4}}{{k^{\prime}}^{4}}y^{*3}{{❘\frac{\partial y}{\partial h}❘}_{\theta_{,{y = y^{*}}}} \cdot {\Delta\theta}}}}},$m = m^(′) = 1or − 1

where k₀ is a radar wavenumber; k and k′ are the wavenumber of two setsof ocean waves which generate the second-order scattering with the radarwave vector; θ is an angle between the ocean wave vector {right arrowover (k)} and the radar wave vector {right arrow over (k₀)}; Δθ is adiscrete interval of θ; S{right arrow over ((⋅))} is a directional wavespectrum; if m=m′=1, σ_(R) ⁽²⁾(ω) indicates the outer second-orderspectrum on the positive first-order peak side; if m=m′=−1, σ_(R) ⁽²⁾(ω)indicates the outer second-order spectrum on the negative first-orderpeak side; Γ is a coupling coefficient; y, y* and h are intermediatevariables defined for the convenience of calculation, wherein y=√{squareroot over (k)} and h=mg^(1/2)y+m′g^(1/2)√{square root over (k′)},wherein g is the acceleration of gravity; and y* is the solution ofconstant Doppler frequency contours denoted by ω−h=0.

(2) obtaining A_(f)(θ) which is defined as the single-frequencycoefficient matrix of the directional wave spectrum by calculating thetheoretical value of a function of the linearized outer second-orderspectrum σ_(RL) ⁽²⁾(ω) divided by the stronger first-order peak energyσ_(R) ⁽¹⁾(ω);

${\frac{\sigma_{RL}^{(2)}(\omega)}{\sigma_{R}^{(1)}\left( \omega_{B} \right)} = {\sum\limits_{\theta = 0}^{2\pi}{{A_{f}(\theta)}{S\left( {m\overset{\rightarrow}{k}} \right)}{\Delta\theta}}}},{{when}{\sigma_{R}^{(1)}(\omega)}{is}{the}{positive}{first} - {order}{peak}}$${\frac{\sigma_{RL}^{(2)}(\omega)}{\sigma_{R}^{(1)}\left( {- \omega_{B}} \right)} = {\sum\limits_{\theta = 0}^{2\pi}{{A_{f}(\theta)}{S\left( {m\overset{\rightarrow}{k}} \right)}{\Delta\theta}}}},{{when}{\sigma_{R}^{(1)}(\omega)}{is}{the}{negative}{first} - {order}{peak}}$

where A_(f)(θ) is the single-frequency coefficient matrix of thedirectional wave spectrum; f in the subscript is the radar frequency;S{right arrow over ((⋅))} is the directional wave spectrum; k is thewavenumber of a set of waves that generate second-order scattering withthe radar wave vector; θ is the angle between the ocean wave vector{right arrow over (k)} and the radar wave vector {right arrow over(k₀)}; Δθ is the discrete interval of θ; σ_(R) ⁽¹⁾(ω_(B)) or σ_(R)⁽¹⁾(−ω_(B)) are the theoretical values of the first-order peaks. Thetheoretical value of the stronger first-order peak is given as:

σ_(R) ⁽¹⁾(ω_(B))=2⁶ πk ₀ ⁴ S(−2{right arrow over (k ₀)}), when σ_(R)⁽¹⁾(ω) is the positive first-order peak

σ_(R) ⁽¹⁾(−ω_(B))=2⁶ πk ₀ ⁴ S(2{right arrow over (k ₀)}), when σ_(R)⁽¹⁾(ω) is the negative first-order peak

whereby the single-frequency coefficient matrix of the directional wavespectrum A_(f)(θ) is given as:

${A_{f}(\theta)} = {{❘\Gamma ❘}^{2}\frac{4\left( {2k_{0}} \right)^{4}y^{*3}}{\left( k^{\prime} \right)^{4}}{❘\frac{\partial y}{\partial h}❘}_{\theta,{y = y^{*}}}}$

where Γ is the coupling coefficient; k₀ is the radar wavenumber; k′ isthe wavenumber of a set of ocean waves which generate the second-orderscattering with the radar wave vector; θ is the angle between the oceanwave vector {right arrow over (k)} and the radar wave vector {rightarrow over (k₀)}; y, y* and h are intermediate variables defined for theconvenience of calculation, wherein y=√{square root over (k)} andh=mg^(1/2)y+m′g^(1/2)√{square root over (k′)}, wherein g is theacceleration of gravity; and y* is the solution of constant Dopplerfrequency contours denoted by ω−h=0.

(3) calculating A_(f) which is defined as the single-frequencycoefficient matrix of the nondirectional wave spectrum by discretizing θinto n pieces at equal intervals and summing all terms of θ in thefollowing equation:

$A_{f} = {\sum\limits_{\theta = 0}^{2\pi}{{A_{f}(\theta)}{G(\theta)}{\Delta\theta}}}$

where A_(f)(θ) is the single-frequency coefficient matrix of thedirectional wave spectrum; f in the subscript is the radar frequency; θis the angle between the ocean wave vector {right arrow over (k)} andthe radar wave vector

$\overset{\rightarrow}{k_{0}};{{\Delta\theta} = \frac{2\pi}{n}}$

is the discrete interval of θ; G(θ) is the direction distributionfunction of the directional wave spectrum.

In a class of this embodiment, the step (f) comprises the step of:

merging R_(f)(ω) which is defined as the function of the outersecond-order spectrum divided by the first-order peak energy for thestronger side of the Doppler spectrum of each frequency and A_(f) whichis defined as the single-frequency coefficient matrix of thenondirectional wave spectrum into a matrix R and a matrix A respectivelyaccording to the following equation:

R=[R _(f) ₁ (ω)R _(f) ₂ (ω)R _(f) ₃ (ω)R _(f) ₄ (ω)]^(T)

A=[A _(f) ₁ A _(f) ₂ A _(f) ₃ A _(f) ₄ ]^(T)

where R_(f) ₁ (ω), R_(f) ₂ (ω), R_(f) ₂ (ω) and R_(f) ₂ (ω) are thefunctions of the outer second-order spectrum to the first-order peakenergy for the stronger side of the Doppler spectrum extracted from thesea echoes of radar frequencies f₁, f₂, f₃, and f₄, respectively; A_(f)₁ , A_(f) ₂ , A_(f) ₃ , and A_(f) ₄ are the coefficient matrices of thenondirectional wave spectrum from the sea echoes of radar frequenciesf₁, f₂, f₃ and f₄, respectively.

In a class of this embodiment, the aforementioned nondirectional wavespectrum in the step (g) is given by the equation:

S(k)=A ⁺ R

Because the matrix A is not square, only the pseudo-inverse of thematrix A can be calculated. A singular value decomposition can becarried out on the matrix A to obtain its pseudo-inverse, and then thenondirectional wave spectrum S(k) can be estimated from the matrix R andthe matrix A⁺.

Compared with the prior art, the method for estimating thenondirectional wave spectrum from the sea echoes of multiple HF radarfrequencies has the superiority that it solves the problem of limitedmeasurable range of wave height observed at a fixed radar frequency bycombining echo signals from multiple radar frequencies, whereby theaccuracy of the method for measuring the nondirectional wave spectrum isimproved and the method can adapt to the complex and changeable seasurface and meet the requirements of robustly extracting wave parametersunder various sea states.

BRIEF DESCRIPTION OF THE DRAWINGS

The disclosure will be more clearly understood from the followingdescription when read in conjunction with the accompanying drawings ofwhich:

FIG. 1 is a flow chart showing the method for directly estimating thenondirectional wave spectrum from the sea echoes of multiple HF radarfrequencies according to an embodiment of the disclosure.

FIGS. 2A-2D are graphs showing the comparison between the nondirectionalwave spectrum estimated from the simulated sea echoes of a single radarfrequency and the theoretical wave spectrum according to an embodimentof the disclosure.

FIG. 3 is a graph showing the comparison between the nondirectional wavespectrum estimated from the simulated sea echoes of multiple radarfrequencies and the theoretical wave spectrum according to an embodimentof the disclosure.

FIGS. 4A-4C are graphs showing the comparison among the nondirectionalwave spectra estimated from the sea echoes collected by amulti-frequency HF radar system according to an embodiment of thedisclosure.

FIG. 5 shows a block diagram of an apparatus for carrying out the methodfor directly estimating the nondirectional wave spectrum from the seaechoes of multiple HF radar frequencies according to an embodiment ofthe disclosure.

DETAILED DESCRIPTION

FIG. 5 shows a block diagram of an apparatus for carrying out a methodfor estimating the nondirectional wave spectrum from the sea echoes ofmultiple HF radar frequencies. The apparatus comprises a multi-frequencyHF radar and a computer processor. The multi-frequency HF radarcomprises a transmitter, a receiver and a signal processor unit. Thetransmitter transmits electromagnetic wave signals of multiplefrequencies between 3-30 MHz in the HF band to the sea surface. Thereceiver receives sea echoes from multiple radar frequencies originatedfrom the electromagnetic wave signals modulated by ocean surfacemovements. The signal processor unit obtains each Doppler spectrum fromthe sea echo of each radar frequency at a fan-shaped unit by performinga first fast Fourier transform (FFT) in distance dimension, a second FFTin Doppler frequency dimension and a digital beamforming. The computerprocessor provides code segments for carrying out the method forestimating the nondirectional wave spectrum. The transmitter andreceiver are respectively connected to the signal processor unit throughcommunication links A and B to realize the full-duplex communication.The signal processor unit is connected to the computer processor throughthe communication link C to realize the full-duplex communication.Communication links A and B include but are not limited to a wirelessconnection and a hardwired connection. Communication link C is a networkconnection, such as Transmission Control Protocol/Internet Protocol(TCP/IP). The multi-frequency HF radar and the computer processor may beinstalled far away from each other.

In operation, the signal processor unit may control the transmitter totransmit frequency modulated interrupted continuous wave (FMICW) chirpsto the ocean surface through the link A. In a Doppler sampling period,the transmitter may transmit FMICW chirps with different frequencies inchronological order, and then the receiver receives sea echoes modulatedby ocean surface movements with corresponding frequencies. The receivercommunicates with the signal processor unit through the link B and sendsthe received sea echo to the signal processor unit. The signal processorunit performs the first FFT in distance dimension and the second FFT inDoppler frequency dimension and the digital beamforming to obtain theDoppler spectrum in a fan-shaped unit. The signal processor unit sendsthe Doppler spectrum in a fan-shaped unit to the computer processorthrough the link C, and then the code segments for estimating thenondirectional wave spectrum provided by the computer processor areexecuted.

FIG. 1 show the outlined steps of the method for estimating thenondirectional wave spectrum from the sea echoes of multiple HF radarfrequencies, comprising the following steps (a)-(g).

The step (a) is dividing a radar detection area into a plurality offan-shaped units at an equal range interval and angle interval accordingto a distance resolution and an angular resolution of a multi-frequencyHF radar, wherein the multi-frequency HF radar is capable ofsimultaneously operating at more than one frequency in the HF band.

The step (b) is obtaining a Doppler spectrum from a sea echo of a singleradar frequency at a fan-shaped unit by performing a first fast FFT indistance dimension, a second FFT in Doppler frequency dimension and adigital beamforming in a signal processor unit, extracting a positivefirst-order peak and a negative first-order peak from the Dopplerspectrum by a peak-searching method and selecting a stronger first-orderpeak σ_(R) ⁽¹⁾(ω) from the positive and negative first-order peaks in acomputer processor.

In the step (b), the Doppler spectrum at the fan-shaped unit is definedas: σ(ω), where ω represents a Doppler frequency generated by the motionof ocean waves to the multi-frequency HF radar; and σ(ω) represents awave energy distribution at the different value of ω.

The positive first-order peak and negative first-order peak areextracted from the Doppler spectrum at the fan-shaped unit using thepeak searching method. The first-order peaks are defined as two peaks inthe Doppler spectrum which are roughly symmetrically distributed on bothsides of zero frequency. The first-order peaks are generated by theBragg scattering of the waves of half the radar wavelength that areeither advancing directly towards the multi-frequency HF radar orreceding directly from the multi-frequency HF radar The detailed stepsof the step (b) are as follows:

searching for a point of a maximum amplitude in the Doppler frequencyrange [0.6ω_(B),1.4ω_(B)] of the Doppler spectrum at the fan-shaped unitas the peak of the positive first-order peak, and recording the Dopplerfrequency of the peak of the positive first-order peak as ω_(P+),wherein ω_(B)=√{square root over (2gk₀)} is a Bragg frequency, wherein

$k_{0} = \frac{2\pi f}{c}$

is a radar wavenumber; c is a speed of light; and f is a radarfrequency;

searching for a local minimum point inside the peak of the positivefirst-order peak where the Doppler frequency meets the inequationω_(P+)−0.2ω_(B)≤ω<ω_(P+), and denoting the Doppler frequency of thelocal minimum point inside the peak of the positive first-order peak asω_(L+);

searching for the local minimum point outside the peak of the positivefirst-order peak where the Doppler frequency satisfies the inequationω_(P+)<(ω) 0.2ω_(B)+ω_(P+), and recording the Doppler frequency of thelocal minimum point outside the peak of the positive first-order peak asω_(R+);

intercepting the Doppler spectrum at the fan-shaped unit with Dopplerfrequency [ω_(L+),ω_(R+)] as the positive first-order peak;

searching for the point of the maximum amplitude in the Dopplerfrequency range [−1.4ω_(B),−0.6ω_(B)] of the Doppler spectrum at thefan-shaped unit as the peak of the negative first-order peak, andrecording the Doppler frequency of the peak of the negative first-orderpeak as ω_(P−);

searching for the local minimum point inside the peak of the negativefirst-order peak where the Doppler frequency meets the inequationω_(P−)<ω≤ω_(P−)+0.2ω_(B), and denoting the Doppler frequency of thelocal minimum point inside the peak of the negative first-order peak asω_(L−);

searching for the local minimum point outside the peak of the negativefirst-order peak where the Doppler frequency satisfies the inequationω_(P−) −0.2ω_(B)≤ω<ω_(P−), and recording the Doppler frequency of thelocal minimum point outside the peak of the negative first-order peak asω_(R−);

intercepting the Doppler spectrum at the fan-shaped unit with Dopplerfrequency [ω_(R−), ω_(L−)] as the negative first-order peak; and

comparing the amplitude of the peak of the positive first-order peakwith that of the negative first-order peak, and selecting thefirst-order peak with a larger amplitude of the peak point as thestronger first-order peak σ_(R) ⁽¹⁾(ω).

The step (c) is dividing a second-order spectrum on the strongerfirst-order peak side into an inner second-order spectrum and an outersecond-order spectrum, and separating the outer second-order spectrum onthe stronger first-order peak side from the Doppler spectrum in thecomputer processor according to a Doppler frequency range of the outersecond-order spectrum;

In the step (c), the second-order spectrum originated from thesecond-order scattering of ocean waves and radar waves is a continuumwith lower amplitude than the first-order peaks and distributed aroundthe first-order peaks.

The step (c) comprises the following sub-steps:

dividing the second-order spectrum on the stronger first-order peak sideinto an inner second-order spectrum and an outer second-order spectrum;

The Doppler frequency range of the outer second-order spectrum on thestronger first-order peak side is given as

$\left\{ \begin{matrix}{{\omega_{c +} < \omega \leq {14\omega_{B}}},{{when}{\sigma_{R}^{(1)}(\omega)}{is}{the}{positive}{first} - {order}{peak}}} \\{{{{- 1.4}\omega_{B}} \leq \omega < \omega_{c -}},{{when}{\sigma_{R}^{(1)}(\omega)}{is}{the}{negative}{first} - {order}{peak}}}\end{matrix} \right.$

The Doppler frequency range of the inner second-order spectrum on thestronger first-order peak side is given as

$\left\{ \begin{matrix}{{{{0.6}\omega_{B}} \leq \omega < \omega_{c +}},{{when}{\sigma_{R}^{(1)}(\omega)}{is}{the}{positive}{first} - {order}{peak}}} \\{{\omega_{c -} < \omega \leq {{- 0.6}\omega_{B}}},{{when}{\sigma_{R}^{(1)}(\omega)}{is}{the}{negative}{first} - {order}{peak}}}\end{matrix} \right.$

where ω is a Doppler frequency generated by the motion of ocean waves tothe multi-frequency HF radar; ω_(B)=√{square root over (2gk₀)} is aBragg frequency; k₀ is the radar wavenumber;

$\omega_{c +} = \frac{\underset{\omega_{L +}}{\sum\limits^{\omega_{R +}}}{\omega \cdot {\sigma_{R}^{(1)}(\omega)}}}{\underset{\omega_{L +}}{\sum\limits^{\omega_{R +}}}{\sigma_{R}^{(1)}(\omega)}}$

is a centroid frequency of the positive first-order peak; and

$\omega_{c -} = \frac{\underset{\omega_{R -}}{\sum\limits^{\omega_{L -}}}{\omega \cdot {\sigma_{R}^{(1)}(\omega)}}}{\underset{\omega_{R -}}{\sum\limits^{\omega_{L -}}}{\sigma_{R}^{(1)}(\omega)}}$

is a centroid frequency of the negative first-order peak.

if the stronger first-order peak obtained in the step (b) is thepositive first-order peak, taking the centroid frequency ω_(c+) as astarting point of Doppler frequency, taking a cutoff frequency of1.4ω_(B) as an endpoint, extracting the Doppler spectrum σ(ω) at thefan-shaped unit in the frequency range of (ω_(c+),1.4ω_(B)], and settingthe Doppler spectrum σ(ω) at the fan-shaped unit in the frequency rangeof (ω_(c+),ω_(R+)] to zero to obtain the outer second-order spectrum onthe stronger first-order peak side σ_(R) ⁽²⁾(ω);

if the stronger first-order peak obtained in the step (b) is thenegative first-order peak, taking the centroid frequency ω_(c−) as thestarting point of Doppler frequency, taking a cutoff frequency of−1.4ω_(B) as the endpoint, extracting the Doppler spectrum σ(ω) at thefan-shaped unit in the frequency range of [−1.4ω_(B),ω_(c−)), andsetting the Doppler spectrum σ(ω) at the fan-shaped unit in thefrequency range of [ω_(R−),ω_(c−)) to zero to obtain the outersecond-order spectrum on the stronger first-order peak side σ_(R)⁽²⁾(ω).

The step (d) is calculating R_(f)(ω) which is a function of the outersecond-order spectrum divided by a first-order peak energy for thestronger side of the Doppler spectrum of the single radar frequency atthe fan-shaped unit in the computer processor. The function R_(f)(ω) canbe given as

$\left\{ \begin{matrix}{{{R_{f}(\omega)} = \frac{\sigma_{R}^{(2)}(\omega)}{\underset{\omega_{L +}}{\sum\limits^{\omega_{R +}}}{{\sigma_{R}^{(1)}(\omega)}{\Delta\omega}}}},{{when}{\sigma_{R}^{(1)}(\omega)}{is}{the}{positive}{first} - {order}{peak}}} \\{{{R_{f}(\omega)} = \frac{\sigma_{R}^{(2)}(\omega)}{\underset{\omega_{R -}}{\sum\limits^{\omega_{L -}}}{{\sigma_{R}^{(1)}(\omega)}{\Delta\omega}}}},{{when}{\sigma_{R}^{(1)}(\omega)}{is}{the}{negative}{first} - {order}{peak}}}\end{matrix} \right.$

where f in a subscript indicates a radar frequency; σ_(R) ⁽¹⁾(ω) is thestronger first-order peak; σ_(R) ⁽²⁾(ω) is the outer second-orderspectrum on the stronger first-order peak side; ω_(L+) is the Dopplerfrequency of the local minimum point inside the peak point of thepositive first-order peak; ω_(R+) is the Doppler frequency of the localminimum point outside the peak point of the positive first-order peak;ω_(L−) is the Doppler frequency of the local minimum point inside thepeak point of the negative first-order peak; φ_(R−) is the Dopplerfrequency of the local minimum point outside the peak point of thenegative first-order peak; and Δω is a frequency resolution of theDoppler spectrum σ(ω).

The step (e) is calculating the single-frequency coefficient matrix ofthe nondirectional wave spectrum by linearizing the outer second-orderspectrum on the side of the stronger first-order peak in the computerprocessor.

The step of calculating the single-frequency coefficient matrix of thenondirectional wave spectrum comprises:

linearizing the outer second-order spectrum on the side of the strongerfirst-order peak;

at the Doppler frequencies near the first-order peak that satisfy thecondition ω_(B)<|ω|≤1.4ω_(B), the direction of one of the two sets ofocean waves that generate the second-order scattering with the radarvector is approximately equal to that of the Bragg wave vector, so theouter second-order spectrum σ_(R) ⁽²⁾(ω) can be linearized as σ_(RL)⁽²⁾(ω)

${{\sigma_{RL}^{(2)}(\omega)} = {\sum\limits_{\theta = 0}^{2\pi}{2^{8}\pi k_{0}^{4}{❘\Gamma ❘}^{2}{S\left( {m\overset{\rightarrow}{k}} \right)}{S\left( {{- 2}m^{\prime}\overset{\rightarrow}{k_{0}}} \right)}\frac{\left( {2k_{0}} \right)^{4}}{k^{\prime 4}}y^{*3}{{❘\frac{\partial y}{\partial h}❘}_{\theta,{y = y^{*}}} \cdot {\Delta\theta}}}}},{m = {m^{\prime} = {{1{or}} - 1}}}$

where

$k_{0} = \frac{2\pi f}{c}$

is a radar wavenumber; k=(0,0.3] and k′=(k²+4kk₀ cosθ+(2k₀)²)^(1/2) arethe wavenumber of two sets of ocean waves which generate thesecond-order scattering with the radar wave vector {right arrow over(k₀)}; θ=[0,2π] is an angle between the ocean wave vector {right arrowover (k)} and the radar wave vector {right arrow over (k₀)}; Δθ is adiscrete interval of θ; S{right arrow over ((⋅))} is a directional wavespectrum; if m=m′=1, σ_(R) ⁽²⁾(ω) indicates the outer second-orderspectrum on the positive first-order peak side; if m=m′=−1, σ_(R) ⁽²⁾(ω)indicates the outer second-order spectrum on the negative first-orderpeak side;

$\Gamma = {{\frac{1}{2}\left\lbrack \frac{{\left( {\overset{\rightarrow}{k} \cdot {\overset{\rightarrow}{k}}_{0}} \right)\left( {{\overset{\rightarrow}{k}}^{\prime} \cdot {\overset{\rightarrow}{k}}_{0}} \right)/k_{0}^{2}} - {2{\overset{\rightarrow}{k} \cdot {\overset{\rightarrow}{k}}^{\prime}}}}{\sqrt{\overset{\rightarrow}{k} \cdot {\overset{\rightarrow}{k}}^{\prime}} - {k_{0}\Delta}} \right\rbrack} - {\frac{i}{2}\left\lbrack {k + k^{\prime} - \frac{\left( {{kk}^{\prime} - {\overset{\rightarrow}{k} \cdot {\overset{\rightarrow}{k}}^{\prime}}} \right)\left( {\omega^{2} + \omega_{B}^{2}} \right)}{{mm}^{\prime}\sqrt{{kk}^{\prime}}\left( {\omega^{2} - \omega_{B}^{2}} \right)}} \right\rbrack}}$

is a coupling coefficient, wherein m, m′, k₀, k, k′, {right arrow over(k)} and {right arrow over (k₀)} are the same as described above, {rightarrow over (k)} is one of two sets of ocean waves which generate thesecond-order scattering with the radar wave vector {right arrow over(k₀)}, ω is the Doppler frequency, ω_(B)=√{square root over (2gk₀)} isthe Bragg frequency, i is an imaginary unit and Δ is the normalizedimpedance of sea surface; y, y* and h are intermediate variables definedfor the convenience of calculation, wherein y=√{square root over (k)}and h=mg_(1/2)y m′g^(1/2)√{square root over (k′)}, wherein g is theacceleration of gravity; and y* is the solution of constant Dopplerfrequency contours denoted by ω−h=0.

(2) obtaining A_(f)(θ) which is defined as the single-frequencycoefficient matrix of the directional wave spectrum by calculating thetheoretical value of a function of the linearized outer second-orderspectrum σ_(R) ⁽²⁾(ω) divided by the stronger first-order peak energyσ_(R) ⁽¹⁾(ω).

${\frac{\sigma_{RL}^{(2)}(\omega)}{\sigma_{R}^{(1)}\left( \omega_{B} \right)} = {\sum\limits_{\theta = 0}^{2\pi}{{A_{f}(\theta)}{S\left( {m\overset{\rightarrow}{k}} \right)}{\Delta\theta}}}},{{when}{\sigma_{R}^{(1)}(\omega)}{is}{the}{positive}{first} - {order}{peak}}$${\frac{\sigma_{RL}^{(2)}(\omega)}{\sigma_{R}^{(1)}\left( {- \omega_{B}} \right)} = {\sum\limits_{\theta = 0}^{2\pi}{{A_{f}(\theta)}{S\left( {m\overset{\rightarrow}{k}} \right)}{\Delta\theta}}}},{{when}{\sigma_{R}^{(1)}(\omega)}{is}{the}{negative}{first} - {order}{peak}}$

where A_(f)(θ) is the single-frequency coefficient matrix of thedirectional wave spectrum; f in the subscript is the radar frequency;S{right arrow over ((⋅))} is the directional wave spectrum; k is thewavenumber of a set of waves that generate second-order scattering withthe radar wave vector; θ=[0,2π] is the angle between the ocean wavevector {right arrow over (k)} and the radar wave vector {right arrowover (k₀)}; Δθ is the discrete interval of θ; σ_(R) ⁽¹⁾(ω_(B)) or σ_(R)⁽¹⁾(−ω_(B)) is the theoretical value of the first-order peaks. Thetheoretical value of the stronger first-order peak is given as:

σ_(R) ⁽¹⁾(ω_(B))=2⁶ πk ₀ ⁴ S(−2{right arrow over (k ₀)}), when σ_(R)⁽¹⁾(ω) is the positive first-order peak

σ_(R) ⁽¹⁾(−ω_(B))=2⁶ πk ₀ ⁴ S(2{right arrow over (k ₀)}), when σ_(R)⁽¹⁾(ω) is the negative first-order peak

whereby the single-frequency coefficient matrix of the directional wavespectrum A_(j)(θ) is given as:

${A_{f}(\theta)} = {{❘\Gamma ❘}^{2}\frac{4\left( {2k_{0}} \right)^{4}y^{*3}}{\left( k^{\prime} \right)^{4}}{❘\frac{\partial y}{\partial h}❘}_{\theta,{y = y^{*}}}}$

where

$\Gamma = {{\frac{1}{2}\left\lbrack \frac{{\left( {\overset{\rightarrow}{k} \cdot {\overset{\rightarrow}{k}}_{0}} \right)\left( {{\overset{\rightarrow}{k}}^{\prime} \cdot {\overset{\rightarrow}{k}}_{0}} \right)/k_{0}^{2}} - {2{\overset{\rightarrow}{k} \cdot {\overset{\rightarrow}{k}}^{\prime}}}}{\sqrt{\overset{\rightarrow}{k} \cdot {\overset{\rightarrow}{k}}^{\prime}} - {k_{0}\Delta}} \right\rbrack} - {\frac{i}{2}\left\lbrack {k + k^{\prime} - \frac{\left( {{kk}^{\prime} - {\overset{\rightarrow}{k} \cdot {\overset{\rightarrow}{k}}^{\prime}}} \right)\left( {\omega^{2} + \omega_{B}^{2}} \right)}{{mm}^{\prime}\sqrt{{kk}^{\prime}}\left( {\omega^{2} - \omega_{B}^{2}} \right)}} \right\rbrack}}$

is the coupling coefficient;

$k_{0} = \frac{2\pi f}{c}$

is the radar wavenumber; k′=(k²+4kk₀ cos θ+(2k₀)²)^(1/2) is thewavenumber of a set of ocean waves which generate the second-orderscattering with the radar wave vector; θ−[0,2π] is the angle between theocean wave vector {right arrow over (k)} and the radar wave vector{right arrow over (k₀)}; y, y* and h are intermediate variables definedfor the convenience of calculation, wherein y=√{square root over (k)}and h=mg^(1/2)y+m′g^(1/2)√{square root over (k′)}, wherein g is theacceleration of gravity; and y* is the solution of constant Dopplerfrequency contours denoted by ω−h=0.

(3) calculating A_(f) which is defined as the single-frequencycoefficient matrix of the nondirectional wave spectrum by discretizing θinto n pieces at equal intervals and summing all terms of θ in thefollowing equation:

$A_{f} = {\sum\limits_{\theta = 0}^{2\pi}{{A_{f}(\theta)}{G(\theta)}{\Delta\theta}}}$

where A_(f)(θ) is the single-frequency coefficient matrix of thedirectional wave spectrum; f in the subscript is the radar frequency;θ−[0,2π] is the angle between the ocean wave vector {right arrow over(k)} and the radar wave vector

${\overset{\rightarrow}{k}}_{0};{{\Delta\theta} = \frac{2\pi}{n}}$

is the discrete interval of θ; G(θ) is the direction distributionfunction of the directional wave spectrum. The cardioid distributionfunction

${G(\theta)} = {{\cos^{2s}\left( \frac{\theta - \theta^{*}}{2} \right)}/{\int_{- \pi}^{\pi}{{\cos^{2s}\left( \frac{\theta}{2} \right)}d\theta}}}$

can be selected, where θ* is the wind direction and s=2 is the spreadfactor.

The step (f) is repeating the above steps (b) to (e); combining the seaechoes of different radar frequencies to construct a matrix R which is afunction of the outer second-order spectrum divided by the first-orderpeak energy for the stronger side of the Doppler spectrum of multipleradar frequencies and merging coefficient matrices of the nondirectionalwave spectrum from multiple radar frequencies into a matrix A in thecomputer processor.

The detailed process of the step (f) is as follows:

R_(f)(ω) which is defined as the function of the outer second-orderspectrum divided by the first-order peak energy for the stronger side ofthe Doppler spectrum of each frequency and A_(f) which is defined as thesingle-frequency coefficient matrix of the nondirectional wave spectrumare merged into a matrix R and a matrix A respectively according to thefollowing equation:

R=[R _(f) ₁ (ω)R _(f) ₂ (ω)R _(f) ₃ (ω)R _(f) ₄ (ω)]^(T)

A=[A _(j) ₁ A _(f) ₂ A _(f) ₃ A _(f) ₄ ]^(T)

where R_(f) ₁ (ω), R_(f) ₂ (ω), R_(f) ₂ (ω) and R_(f) ₂ (ω) are thefunctions of the outer second-order spectrum divided by the first-orderpeak energy for the stronger side of the Doppler spectrum extracted fromthe sea echoes of radar frequencies f₁, f₂, f₃ and f₄, respectively; andA_(f) ₁ , A_(f) ₂ , A_(f) ₃ , and A_(f) ₄ are the coefficient matricesof the nondirectional wave spectrum from the sea echoes of radarfrequencies fj, f₃ and f₄ respectively.

The step (g) is calculating the pseudo-inverse A⁺ of the matrix A by asingular value decomposition; and estimating the nondirectional wavespectrum at the fan-shaped unit from the matrix R and the pseudo-inverseA⁺ in the computer processor.

The aforementioned nondirectional wave spectrum in the step (g) is givenby the equation:

S(k)=A ⁺ R

Because the matrix A is not square, only the pseudo-inverse of thematrix A can be calculated. The singular value decomposition can becarried out on the matrix A to obtain its pseudo-inverse, and then thenondirectional wave spectrum S (k) can be estimated from the matrix Rand the matrix A⁺.

As an example, the Pierson-Moskowitz (PM) spectrum and the cardioiddistribution function are selected to simulate the Doppler spectrum ofthe sea echoes of radar frequency 8 MHz, 13 MHz, 19 MHz and 25 MHzrespectively. The wind speed is set to 10 m/s, and the wind direction isset to 135°. The comparison between the nondirectional wave spectrumestimated from the simulated sea echoes of a single radar frequency andthe theoretical wave spectrum is shown in FIGS. 2A-2D. As shown in FIG.3 the nondirectional wave spectrum estimated from the sea echoes ofmultiple radar frequencies is in good agreement with that of thetheoretical values. FIGS. 2A-2D and FIG. 3 demonstrate that theestimation result of the nondirectional wave spectrum from the multiplefrequencies provides reasonable agreement with the theoretical valuescompared with that from the sea echoes at 8 MHz, 13 MHz, 19 MHz, and 25MHz.

The method provided by the disclosure is applied to a dataset collectedby a multi-frequency HF radar system operating at 8.267 MHz and 19.2MHz. The wave height estimated by the method from the radar data iscompared with that measured by a buoy to evaluate the performance of themethod. The comparison among the nondirectional wave spectra estimatedfrom the sea echoes collected with a multi-frequency HF radar system at08:00 AM on Jul. 10, 2015 is shown in FIGS. 4A-4C. The wave heightextracted by the nondirectional wave spectrum directly estimated fromthe multifrequency data is more consistent with that measured by thebuoy than that estimated from the single-frequency data.

The method provided by the disclosure for directly estimating thenondirectional wave spectrum from the sea echoes of multiple HF radarfrequencies breaks through the limitation of the sea states that can bedetected at a fixed radar frequency by combining echo signals frommultiple radar frequencies, whereby the robustness and accuracy ofestimating the nondirectional wave spectrum are greatly enhanced and thepotential of the method for monitoring the complicated and changeablesea surface is recognized. Therefore, it is feasible to estimate thenondirectional wave spectrum from the sea echoes of multiple HF radarfrequencies, and it has higher accuracy and wide application prospectcompared with the conventional method for inverting the nondirectionalwave spectrum from the single-frequency sea echo.

It will be obvious to those skilled in the art that changes andmodifications may be made, and therefore, the aim in the appended claimsis to cover all such changes and modifications.

What is claimed is:
 1. A method for directly estimating a nondirectional wave spectrum from sea echoes of multiple HF radar frequencies, the method comprising: (a) dividing a radar detection area into a plurality of fan-shaped units at an equal range interval and angle interval according to a distance resolution and an angular resolution of a multi-frequency HF radar, wherein the multi-frequency HF radar is capable of simultaneously operating at more than one frequency in the HF band; (b) obtaining a Doppler spectrum from a sea echo of a single radar frequency at a fan-shaped unit by performing a first fast Fourier transform (FFT) in distance dimension, a second FFT in Doppler frequency dimension and a digital beamforming in a signal processor unit, extracting a positive first-order peak and a negative first-order peak from the Doppler spectrum by a peak-searching method and selecting a stronger first-order peak σ_(R) ⁽¹⁾(ω) from the positive and negative first-order peaks in a computer processor; (c) dividing a second-order spectrum on a stronger first-order peak side into an inner second-order spectrum and an outer second-order spectrum and separating the outer second-order spectrum on the stronger first-order peak side from the Doppler spectrum in the computer processor according to a Doppler frequency range of the outer second-order spectrum; (d) calculating R_(f)(ω) which is a function of the outer second-order spectrum divided by a first-order peak energy for the stronger side of the Doppler spectrum of the single radar frequency at the fan-shaped unit in the computer processor; (e) calculating A_(f) which is defined as a single-frequency coefficient matrix of the nondirectional wave spectrum by linearizing the outer second-order spectrum on the side of the stronger first-order peak in the computer processor; (f) repeating the above steps (b) to (e), combining the sea echoes of different radar frequencies to construct a matrix R which is a function of the outer second-order spectrum divided by the first-order peak energy for the stronger side of the Doppler spectrum of multiple radar frequencies and merging coefficient matrices of the nondirectional wave spectrum from multiple radar frequencies into a matrix A in the computer processor; and (g) calculating a pseudo-inverse A⁺ of the matrix A by a singular value decomposition and estimating the nondirectional wave spectrum at the fan-shaped unit from the matrix R and the pseudo-inverse A⁺ in the computer processor.
 2. The method of claim 1, wherein in (b), the Doppler spectrum at the fan-shaped unit is defined as: σ(ω), where ω represents a Doppler frequency generated by the motion of ocean waves to the multi-frequency HF radar; and σ(ω) represents a wave energy distribution at a different value of ω; wherein in (b), the positive first-order peak and negative first-order peak are extracted from the Doppler spectrum at the fan-shaped unit using the peak searching method; the first-order peaks are defined as two peaks in the Doppler spectrum which are roughly symmetrically distributed on both sides of zero frequency; the first-order peaks are generated by the Bragg scattering of the waves of half the radar wavelength that are either advancing directly towards the multi-frequency HF radar or receding directly from the multi-frequency HF radar; wherein (b) comprises: searching for a point of a maximum amplitude in the Doppler frequency range [0.6ω_(B),1.4ω_(B)] of the Doppler spectrum at the fan-shaped unit as the peak of the positive first-order peak, and recording the Doppler frequency of the peak of the positive first-order peak as ω_(P+), wherein ω_(B)=√{square root over (2gk₀)}, is a Bragg frequency, wherein $k_{0} = \frac{2\pi f}{c}$ is a radar wavenumber; c is a speed of light; and f is a radar frequency; searching for a local minimum point inside the peak of the positive first-order peak where the Doppler frequency meets the inequation ω_(P+)−0.2ω_(B)≤ω<ω_(P+), and denoting the Doppler frequency of the local minimum point inside the peak of the positive first-order peak as ω_(L+); searching for the local minimum point outside the peak of the positive first-order peak where the Doppler frequency satisfies the inequation ω_(P+)<ω≤0.2ω_(B)+ω_(P+), and recording the Doppler frequency of the local minimum point outside the peak of the positive first-order peak as ω_(R+); intercepting the Doppler spectrum at the fan-shaped unit with Doppler frequency [ω_(L+), ω_(R+)] as the positive first-order peak; searching for the point of the maximum amplitude in the Doppler frequency range [−1.4ω_(B), −0.6ω_(B)] of the Doppler spectrum at the fan-shaped unit as the peak of the negative first-order peak, and recording the Doppler frequency of the peak of the negative first-order peak as ω_(P−); searching for the local minimum point inside the peak of the negative first-order peak where the Doppler frequency meets the inequation ω_(P−)<ω≤ω_(P−)+0.2ω_(B), and denoting the Doppler frequency of the local minimum point inside the peak of the negative first-order peak as ω_(L−); searching for the local minimum point outside the peak of the negative first-order peak where the Doppler frequency satisfies the inequation ω_(P−)−0.2φ_(B)≤ω<ω_(P−), and recording the Doppler frequency of the local minimum point outside the peak of the negative first-order peak as ω_(R−); intercepting the Doppler spectrum at the fan-shaped unit with Doppler frequency [ω_(R−), ω_(L−)] as the negative first-order peak; and comparing the amplitude of the peak of the positive first-order peak with that of the negative first-order peak, and selecting the first-order peak with a larger amplitude of the peak point as the stronger first-order peak σ_(R) ⁽¹⁾(ω).
 3. The method of claim 1, wherein in (c), the second-order spectrum originated from the second-order scattering of ocean waves and radar waves is a continuum with lower amplitude than the first-order peaks and distributed around the first-order peaks; wherein (c) comprises the steps of: (1) dividing the second-order spectrum on the stronger first-order peak side into an inner second-order spectrum and an outer second-order spectrum; the Doppler frequency range of the outer second-order spectrum on the stronger first-order peak side is given as $\left\{ {\begin{matrix} {{\omega_{c +} < \omega \leq {1.4\omega_{B}}}\ ,{{when}{\sigma_{R}^{(1)}(\omega)}{is}{the}{positive}{first} - {order}{peak}}} \\ {{{{- 1.4}\omega_{B}} \leq \omega < \omega_{c -}},{{when}{\sigma_{R}^{(1)}(\omega)}{is}{the}{negative}{first} - {order}{peak}}} \end{matrix};} \right.$ the Doppler frequency range of the inner second-order spectrum on the stronger first-order peak side is given as $\left\{ {\begin{matrix} {{{0.6\omega_{B}} \leq \omega < \omega_{c +}},{{when}{\sigma_{R}^{(1)}(\omega)}{is}{the}{positive}{first} - {order}{peak}}} \\ {{\omega_{c -} < \omega \leq {{- 0.6}\omega_{B}}},{{when}{\sigma_{R}^{(1)}(\omega)}{is}{the}{negative}{first} - {order}{peak}}} \end{matrix};} \right.$ where ω is a Doppler frequency generated by the motion of ocean waves to a radar; ω_(B)=√{square root over (2gk₀)} is a Bragg frequency; k₀ is a radar wavenumber; ${\omega_{c +} = {\frac{\overset{\omega_{R +}}{\sum\limits_{\omega_{L +}}}{\omega \cdot {\sigma_{R}^{(1)}(\omega)}}}{\overset{\omega_{R +}}{\sum\limits_{\omega_{L +}}}{\sigma_{R}^{(1)}(\omega)}}{is}a{centroid}{frequency}{of}{the}{positive}{first} - {order}{peak}}};{and}$ ${\omega_{c -} = {\frac{\overset{\omega_{L -}}{\sum\limits_{\omega_{R -}}}{\omega \cdot {\sigma_{R}^{(1)}(\omega)}}}{\overset{\omega_{L -}}{\sum\limits_{\omega_{R -}}}{\sigma_{R}^{(1)}(\omega)}}{is}a{centroid}{frequency}{of}{the}{negative}{first} - {order}{peak}}};$ (2) if the stronger first-order peak obtained in (b) of claim 1 is the positive first-order peak, taking the centroid frequency ω_(c+), as a starting point of Doppler frequency, taking a cutoff frequency of 1.4ω_(B) as an endpoint, extracting the Doppler spectrum at the fan-shaped unit σ(ω) in the frequency range of (ω_(c+),1.4ω_(B)], and setting the Doppler spectrum at the fan-shaped unit σ(ω) in the frequency range of (ω_(c+),ω_(R+)] to zero to obtain the outer second-order spectrum on the stronger first-order peak side σ_(R) ⁽²⁾(ω); if the stronger first-order peak obtained in (b) of claim 1 is the negative first-order peak, taking the centroid frequency ω_(c−) as the starting point of Doppler frequency, taking a cutoff frequency of −1.4ω_(B) as the endpoint, extracting the Doppler spectrum at the fan-shaped unit σ(ω) in the frequency range of [−1.4ω_(B),ω_(c−)), and setting the Doppler spectrum at the fan-shaped unit σ(ω) in the frequency range of [ω_(R−),ω_(c−)) to zero to obtain the outer second-order spectrum on the stronger first-order peak side σ_(R) ⁽²⁾(ω).
 4. The method of claim 1, wherein in (d), the aforementioned function R_(f)(ω) is given by the equation: $\left\{ {\begin{matrix} {{{R_{f}(\omega)} = \frac{\sigma_{R}^{(2)}(\omega)}{\overset{\omega_{R +}}{\sum\limits_{\omega_{L +}}}{{\sigma_{R}^{(1)}(\omega)}{\Delta\omega}}}},{{when}{\sigma_{R}^{(1)}(\omega)}{is}{the}{positive}{first} - {order}{peak}}} \\ {{{R_{f}(\omega)} = \frac{\sigma_{R}^{(2)}(\omega)}{\overset{\omega_{L -}}{\sum\limits_{\omega_{R -}}}{{\sigma_{R}^{(1)}(\omega)}{\Delta\omega}}}},{{when}{\sigma_{R}^{(1)}(\omega)}{is}{the}{negative}{first} - {order}{peak}}} \end{matrix};} \right.$ where f in a subscript indicates a radar frequency; σ_(R) ⁽¹⁾(ω) is the stronger first-order peak; σ_(R) ⁽²⁾ (ω) is the outer second-order spectrum on the stronger first-order peak side; ω_(L+) is the Doppler frequency of the local minimum point inside the peak of the positive first-order peak; ω_(R+) is the Doppler frequency of the local minimum point outside the peak of the positive first-order peak; ω_(L−) is the Doppler frequency of the local minimum point inside the peak of the negative first-order peak; ω_(R−) is the Doppler frequency of the local minimum point outside the negative first-order peak; and Δω is a frequency resolution of the Doppler spectrum σ(ω).
 5. The method of claim 1, wherein (e) comprises the steps of: (1) linearizing the outer second-order spectrum on the side of the stronger first-order peak; at the Doppler frequencies near the first-order peak that satisfy the condition ω_(B)<|ω|≤1.4ω_(B), the direction of one of the two sets of ocean waves that generate the second-order scattering with the radar vector is approximately equal to that of the Bragg wave vector, so the outer second-order spectrum σ_(R) ⁽²⁾ (ω) is linearized as σ_(RL) ⁽²⁾(ω): ${{\sigma_{RL}^{(2)}(\omega)} = {\sum\limits_{\theta = 0}^{2\pi}{2^{8}\pi k_{0}^{4}{❘\Gamma ❘}^{2}{S\left( {m\overset{\rightarrow}{k}} \right)}{S\left( {{- 2}m^{\prime}\overset{\rightarrow}{k_{0}}} \right)}\frac{\left( {2k_{0}} \right)^{4}}{k^{\prime 4}}y^{*3}{{❘\frac{\partial y}{\partial h}❘}_{\theta,{y = y^{*}}} \cdot {\Delta\theta}}}}},{m = {m^{\prime} = {{1{or}} - 1}}}$ where k₀ is a radar wavenumber; k and k′ are the wavenumber of two sets of ocean waves which generate the second-order scattering with the radar wave vector {right arrow over (k₀)}; θ is an angle between the ocean wave vector {right arrow over (k)} and the radar wave vector {right arrow over (k₀)}; Δθ is a discrete interval of θ; S{right arrow over ((⋅))} is a directional wave spectrum; if m=m′=1, σ_(R) ⁽²⁾(ω) indicates the outer second-order spectrum on the positive first-order peak side; if m=m′=−1, σ_(R) ⁽²⁾(ω) indicates the outer second-order spectrum on the negative first-order peak side; ω is a Doppler frequency; Γ is a coupling coefficient; y, y* and h are intermediate variables defined for the convenience of calculation, wherein y=√{square root over (k)} and h+mg^(1/2)y+m′g^(1/2)√{square root over (k′)}, g is the acceleration of gravity; and y* is the solution of constant Doppler frequency contours denoted by ω−h=0; (2) obtaining A_(f) (θ) which is defined as the single-frequency coefficient matrix of the directional wave spectrum by calculating the theoretical value of a function of the linearized outer second-order spectrum σ_(RL) ⁽²⁾(ω) divided by the stronger first-order peak energy σ_(R) ⁽¹⁾(ω); ${\frac{\sigma_{RL}^{(2)}(\omega)}{\sigma_{R}^{(1)}\left( \omega_{B} \right)} = {\sum\limits_{\theta = 0}^{2\pi}{{A_{f}(\theta)}{S\left( {m\overset{\rightarrow}{k}} \right)}{\Delta\theta}}}},{{when}{\sigma_{R}^{(1)}(\omega)}{is}{the}{positive}{first} - {order}{peak}}$ ${\frac{\sigma_{RL}^{(2)}(\omega)}{\sigma_{R}^{(1)}\left( {- \omega_{B}} \right)} = {\sum\limits_{\theta = 0}^{2\pi}{{A_{f}(\theta)}{S\left( {m\overset{\rightarrow}{k}} \right)}{\Delta\theta}}}},{{when}{\sigma_{R}^{(1)}(\omega)}{is}{the}{negative}{first} - {order}{peak}}$ where A_(f) (θ) is the single-frequency coefficient matrix of the directional wave spectrum; f in the subscript is the radar frequency; S{right arrow over ((⋅))} is the directional wave spectrum; k is the wavenumber of a set of waves that generate second-order scattering with the radar wave vector; θ is the angle between the ocean wave vector {right arrow over (k)} and the radar wave vector {right arrow over (k₀)}; Δθ is the discrete interval of θ; σ_(R) ⁽¹⁾(ω_(B)) and σ_(R) ⁽¹⁾(−ω_(B)) are the theoretical values of the first-order peaks; the theoretical value of the stronger first-order peak is given as: σ_(R) ⁽¹⁾(ω_(B))=2⁶ πk ₀ ⁴ S(−2{right arrow over (k ₀)}), when σ_(R) ⁽¹⁾(ω) is the positive first-order peak σ_(R) ⁽¹⁾(−ω_(B))=2⁶ πk ₀ ⁴ S(2{right arrow over (k ₀)}), when σ_(R) ⁽¹⁾(ω) is the negative first-order peak whereby the single-frequency coefficient matrix of the directional wave spectrum A_(j)(θ) is given as: ${A_{f}(\theta)} = {{❘\Gamma ❘}^{2}\frac{4\left( {2k_{0}} \right)^{4}y^{*3}}{\left( k^{\prime} \right)^{4}}{❘\frac{\partial y}{\partial h}❘}_{\theta,{y = y^{*}}}}$ where Γ is the coupling coefficient; k₀ is the radar wavenumber; k′ is the wavenumber of a set of ocean waves which generate the second-order scattering with the radar wave vector; θ is the angle between the ocean wave vector {right arrow over (k)} and the radar wave vector {right arrow over (k₀)}; y, y* and h are intermediate variables defined for the convenience of calculation, wherein y=√{square root over (k)} and h=mg^(1/2)y+m′g^(1/2){right arrow over (k′)}, wherein g is the acceleration of gravity; and y* is the solution of constant Doppler frequency contours denoted by ω−h=0; (3) calculating A_(f) which is defined as the single-frequency coefficient matrix of the nondirectional wave spectrum by discretizing θ into n pieces at equal intervals and summing all terms of θ in the following equation: $A_{f} = {\sum\limits_{\theta = 0}^{2\pi}{{A_{f}(\theta)}{G(\theta)}{\Delta\theta}}}$ where A_(f)(θ) is the single-frequency coefficient matrix of the directional wave spectrum; f in the subscript is the radar frequency; θ is the angle between the ocean wave vector {right arrow over (k)} and the radar wave vector ${\overset{\rightarrow}{k}}_{0};{{\Delta\theta} = \frac{2\pi}{n}}$ is the discrete interval of θ; G(θ) is the direction distribution function of the directional wave spectrum.
 6. The method of claim 1, wherein (f) comprises the step of: merging R_(f)(ω) which is defined as the function of the outer second-order spectrum divided by the first-order peak energy for the stronger side of the Doppler spectrum of each frequency and A_(f) which is defined as the single-frequency coefficient matrix of the nondirectional wave spectrum into a matrix R and a matrix A respectively according to the following equation: R=[R _(f) ₁ (ω)R _(f) ₂ (ω)R _(f) ₃ (ω)R _(f) ₄ (ω)]^(T) A=[A _(f) ₁ A _(f) ₂ A _(f) ₃ A _(f) ₄ ]^(T) where R_(f) ₁ (ω), R_(f) ₂ (ω), R_(f) ₂ (ω) and R_(f) ₂ (ω) are the functions of the outer second-order spectrum divided by the first-order peak energy for the stronger side of the Doppler spectrum extracted from the sea echoes of radar frequencies f₁, f₂, f₃ and f₄, respectively; A_(f) ₁ , A_(f) ₂ , A_(f) ₃ , and A_(f) ₄ are the coefficient matrices of the nondirectional wave spectrum from the sea echoes of radar frequencies f₁, f₂, f₃ and f₄, respectively.
 7. The method of claim 1, wherein (g) comprises the steps of: calculating the pseudo-inverse A⁺ of the matrix A by the singular value decomposition; and estimating the nondirectional wave spectrum S(k) from the matrix R and the pseudo-inverse A⁺ by the equation S(k)=A⁺R.
 8. The method of claim 1, being carried out in an apparatus, the apparatus comprising a multi-frequency HF radar and a computer processor: the multi-frequency HF radar comprising: a transmitter for transmitting electromagnetic wave signals with multiple frequencies between 3-30 MHz in the HF band to the sea surface; a receiver for receiving sea echoes of multiple radar frequencies originated from the electromagnetic wave signals modulated by sea surface movements; and a signal processor unit for obtaining each Doppler spectrum from the sea echo of each radar frequency at a fan-shaped unit by performing the first FFT in distance dimension, the second FFT in Doppler frequency dimension and the digital beamforming; the computer processor providing code segments for: (1) extracting a positive first-order peak and a negative first-order peak from a Doppler spectrum of a single radar frequency by a peak-searching method and selecting a stronger first-order peak σ_(R) ⁽¹⁾(ω) from the positive and negative first-order peaks; (2) dividing a second-order spectrum on a stronger first-order peak side into an inner second-order spectrum and an outer second-order spectrum and separating the outer second-order spectrum on the stronger first-order peak side from the Doppler spectrum according to a Doppler frequency range of the outer second-order spectrum; (3) calculating R_(f)(ω) which is a function of the outer second-order spectrum divided by a first-order peak energy for the stronger side of the Doppler spectrum of the single radar frequency at the fan-shaped unit; (4) calculating A_(f) which is defined as a single-frequency coefficient matrix of the nondirectional wave spectrum by linearizing the outer second-order spectrum on the side of the stronger first-order peak; (5) repeating the above steps (1) to (4), combining the sea echoes of different radar frequencies to construct a matrix R which is a function of the outer second-order spectrum divided by the first-order peak energy for the stronger side of the Doppler spectrum of multiple radar frequencies and merging coefficient matrices of the nondirectional wave spectrum from multiple radar frequencies into a matrix A; and (6) calculating a pseudo-inverse A⁺ of the matrix A by a singular value decomposition and estimating the nondirectional wave spectrum at the fan-shaped unit from the matrix R and the pseudo-inverse A⁺. 